#########################################################################################################
# Part 5: Prior Density, Posterior Density and Relative Belief Ratio for the Parameter of Interest psi  #
#########################################################################################################
# May 1, 2024


# psi is the quantity we want to make inference about and is defined as a function of mu, Sigma and xi
psifn = function(muval, xival, Sigmaval){
# muval is a vector of means
# xival is a precision matrix associated with variance matrix Sigmaval
# the code here depends on the psi function we wish to make inference about
      psi=muval[1]
      return(psi)
}


########################################################################################################
# NOTE: to use this program successfully there needs to be some iteration as follows                   #
#1. generate the prior sample of Nprior using Part 2 and here obtain min and max sample values of psi  #
#2. in this program choose numcells = number of cells in a density histogram for both prior            # 
#   and posterior and divide the interval (min, max) into numcells equql length subintervals           #
#3. plot the prior density and smooth as needed although if too much smoothing is required it is       #
#   probably better to increase Nprior or decrease numcells                                            #
#4. compute the importance sampled values of psi                                                       #
#5. check that the importance sampled values are covered by (min, max), if not there is probably       # 
#   prior-data conflict                                                                                #
#6. get the density histogram estimate of the posterior density of psi using the importance            # 
#   sampled values of psi but note that this is estimated differently than for the prior               #
#7. if the estimate of the posterior density is too rough and smoothing makes it too degenerate        #
#   then the solution is to increase numcells and possibly both Nprior and Npostimp                    #
########################################################################################################

##########################################################################################################
# obtain the prior density of psi

# obtain prior sample generated in Part 2
mu_prior = sample_prior_vals$mu_matrix
Sigma_prior =  sample_prior_vals$Sigma_mat[]
xi_prior = sample_prior_vals$xi_mat


# compute prior sample of psi values 
prior_psi_vals = numeric(Nprior)
for (i in 1:Nprior){
    prior_psi_vals[i] = psifn(mu_prior[i],xi_prior[[i]],Sigma[[i]])
}

# display a density histogram of the values of prior generated from its prior and see what the range of values is
hist(prior_psi_vals,freq=F)
prior_psi_lower_bd = min(prior_psi_vals)
prior_psi_upper_bd = max(prior_psi_vals)
prior_psi_lower_bd
prior_psi_upper_bd

# now break up the range into numcells subintervals of equal length and plot the density histogram
numcells = 100 # this has to be input by the user
delta_psi = (prior_psi_upper_bd - prior_psi_lower_bd) /numcells 
breaks = seq(prior_psi_lower_bd, prior_psi_upper_bd, by = delta_psi)
delta_psi
breaks
hist(prior_psi_vals, breaks, freq = F)

# if the previous density histogram looks reasonable then get the density values and mid points of the cells
# otherwise modify numbreaks
prior_psi_hist = hist(prior_psi_vals, breaks, freq = F)
prior_psi_mids = prior_psi_hist$mids
prior_psi_density = prior_psi_hist$density 

# smoothed plot of the prior density of psi
# if a plot is too rough smooth by averaging
# mprior = an odd number of points to average prior density values, (mprior=1,3,5,...) by averaging
# the density value, (mprior-1)/2 values to the left and (mprior-1)/2 values to the right
prior_psi_dens_smoothed = prior_psi_density
mprior=7
halfm=(mprior-1)/2
for (i in (1+halfm):(numcells-halfm)){
    sum = 0
    for (j in (-halfm):halfm){
    sum = sum + prior_psi_density[i+j]
    }
prior_psi_dens_smoothed[i]=sum/mprior  
}
plot(prior_psi_mids, prior_psi_dens_smoothed, type="l", xlab="psi", ylab="density", main="The prior density of psi")

#####################################################################################################
# obtain posterior density of psi

# obtain importance sample generated in Part 3
imp_mu = imp_vals$mu_xi
imp_Sigma=imp_vals$Sigma
imp_xi = imp_vals$xi
imp_weights = imp_vals$weights_vector

# compute sample of psi values generated by the importance sampler 
imp_psi_vals = numeric(Npostimp)
for (i in 1:Npostimp){
    imp_psi_vals[i] = psifn(imp_mu[i],imp_xi[[i]],imp_Sigma[[i]])
}

# compute the estimate of the posterior cdf of psi
nbreaks=numcells+1
post_psi_cdf = rep(0, nbreaks)
for(i in 2:nbreaks){
  for(j in 1:Npostimp){
    if(imp_psi_vals[j] <= breaks[i]){
      post_psi_cdf[i] = post_psi_cdf[i] + imp_weights[j]
    }
  }
}

# check interval (prior_psi_lower_bd, prior_psi_upper_bd) covers importance sampling values of psi.
min(imp_psi_vals) - breaks[1] # should be positive
max(imp_psi_vals) - breaks[nbreaks] # should be negative

# compute posterior density of psi
post_psi_density = diff(post_psi_cdf)/delta_psi

# smoothed plot of the posterior density of psi
# if a plot is too rough smooth by averaging
# mprior = an odd number of points to average prior density values, (mprior=1,3,5,...) by averaging
# the density value, (mpost-1)/2 values to the left and (mpost-1)/2 values to the right
post_psi_dens_smoothed = post_psi_density
mpost = 5
halfm = (mpost-1)/2
for (i in (1+halfm):(numcells-halfm)){
    sum = 0
    for (j in (-halfm):halfm){
    sum = sum + post_psi_density[i+j]
    }
post_psi_dens_smoothed[i] = sum/mpost  
}
plot(prior_psi_mids, post_psi_dens_smoothed, type="l", xlab="psi", ylab="density", main="The posterior density of psi")

####################################################################################################
# obtain the relative belief ratio of psi
RB_psi = rep(0,numcells)
for (i in 1:numcells){
if (prior_psi_dens_smoothed[i] != 0){
RB_psi[i] = post_psi_dens_smoothed[i]/prior_psi_dens_smoothed[i]}
}
plot(prior_psi_mids, RB_psi, type="l", xlab="psi", ylab="RBR", main="The relative belief ratio of psi")

#############################################################################################################
# The plots side-by-side
par(mfrow = (c(1, 3)))
plot(prior_psi_mids, prior_psi_dens_smoothed, type = "l", xlab = "psi", ylab = "density", col = "red", main="The prior density of psi")
plot(prior_psi_mids, post_psi_dens_smoothed, type = "l", xlab = "psi", ylab = "density", col = "blue", main="The posterior density of psi" )
plot(prior_psi_mids, RB_psi, type = "l", xlab = "psi", ylab = "RBR", col = "green", main="The relative belief ratio of psi")